Mixed Finite Elements for the Nonlinear Richards’ Equation

In this work, we present a two-dimensional mixed-hybrid finite element model of variably saturated flow on unstructured triangular meshes. Velocities are approximated using lowest order Raviart-Thomas (RT0) elements with piecewise constant pressure. The resulting nonlinear systems of algebraic equations are solved using Picard or Newton iterations in combination with ad hoc preconditioning techniques to improve the convergence of the linearized mixed system which is known to be ill-conditioned particularly in the case of stationary flow. Two sample tests show that the Newton approach achieves fast convergence in many situations if a good initial guess is provided by either the Picard scheme or relaxation techniques. 1 The Mathematical Model Richard’s equation is commonly used in modeling the flow of water in the unsaturated zone. It is obtained by combining Darcy’s law with the continuity equation. We write Richard’s equation using pressure head = p= , where p is pressure and is the specific weight of the fluid, as primary variable [1]. The water flux ~q in variably saturated porous media is governed by the extended form of Darcy’s law [2, 3]: ~q = Kr( )Ks~ r ( + z) (1) where Kr( ) is the relative hydraulic conductivity, Ks is the saturated hydraulic conductivity tensor, ~ r is the gradient operator, and z is the vertical coordinate directed upward. The mass conservation equation for variably saturated flow can be written as: (Sw)@ @t + ~ r ~q = f (2) where (Sw) = SwSs+n@Sw=@ is the general storage term, Sw is the water saturation, Ss the elastic storage coefficient, n is the porosity of the medium, t is time, and f represents the source or sink term. Equation (1) is nonlinear because of the dependency of the relative conductivity on pressure head, while (2) assumes its nonlinearity in the general storage term, a function of pressure head through the Sw relationship. The nonlinear functions Kr( ) and Sw( ) in (1) and (2), respectively, can be modeled using different characteristic relationships that are generally determined experimentally for different soils. Neuman and Dirichlet boundary conditions complete the mathematical formulation. Under certain conditions, for example large infiltration on very dry soils, equation (2) becomes highly nonlinear, and sharp gradients appear in its solution. In this case, numerical schemes may show extremely slow convergence rates, or even non convergence [4], unless very small time steps are employed. 2 Mixed-hybrid discretization Let T be a triangulation of the domain R 2 with m triangles (Tj; j = 1 : : :m) and n edges (ej; j = 1; : : : n). The pressure head and the water flux ~q can be approximated by: = m X j=1 pj j ~q = 3m X j=1 qj ~ wj where j and ~ wj are scalar and vector basis functions. We use the Raviart-Thomas space of degree zero (RT0) [5] whereby the functions i are P0 polynomials equal to one on triangle and zero elsewhere, while ~ wj are P1 polynomials in each coordinate direction, defined on the edge ej of triangle Ti by the condition: Zei ~ wj ~nids = ij i; j = 1; : : : ; 3m The support of ~ wj is given by one of the two triangles sharing edge ej . The mixed-hybrid formulation defines the flux ~q as a function with discontinuous normal components across the inter-element boundaries. The desired continuity is imposed on ~q through the definition of one extra set of unknowns, the so called Lagrange multipliers [6]: = n0 Xj j j where n0 is the number of internal edges and the basis function j is equal to 1 over ej and 0 elsewhere. The hybrid formulation can be written as: Aq B + C = g BTq + D@ @t = f CTq = 0 where A (aij) = m X̀=1 ZT` K` 1 ~ wì ~ wj̀d B (bij) = m X̀=1 ZT` ~ r ~ wì j̀d C (cij) = m X̀=1 I@T` j̀ ~ wì ~nì d D (dij) = m X̀=1 ZT` ` ì j̀d Vectors g and f take into account boundary and source terms. For the time discretization we use a weighted finite difference scheme with 2]0; 1] as the weighting parameter and define 0B@ q 1CAk+ 0B@ q 1CAk+1 + (1 )0B@ q 1CAk With the notation Ak+ and Dk+ we mean that A and D, the only matrices which vary depending on pressure head, are evaluated at k+ . The final system of nonlinear equations can be written as: ( E +G)k+ uk+1 = ( (1 )E +G)k+ uk + r (3) E = 0B@ A B C BT 0 0 CT 0 0 1CA ; G = 0B@ 0 0 0 0 1 tD 0 0 0 0 1CA ; uk+1 = 0B@ qk+1 k+1 k+1 1CA ; r = 0B@ g f0 1CA 3 Newton-like linearizations The solution of the nonlinear system (3) is obtained by means of Newton-like iterative schemes [7]. Define the the general nonlinear function of uk+1 as: F (uk+1) = 0 where F (uk+1) = Ek+ uk+ +Gk+ (uk+1 uk) r The Newton iteration can be generally defined as: J(um)sm = F (um) um+1 = um + sm where J = (@Fi=@uk+1 j ) is the Jacobian matrix, m is the iteration index, the vector sm = um+1 um is the Newton, or search, direction, and um um;k+1 and um+1 um+1;k+1. The Jacobian matrix can be written explicitly as: J = @F @u = ( E +G) + @E @uk+1uk+ + @G @uk+1 uk+1 uk Rearranging the terms we obtain: ( E +G+ T ) sm = E + G um; + G uk + r (4) um+1 = um + sm where T = @E @uum; + @G @u um uk Here and in the following, all the matrices are considered as calculated at um; = um;k+ . An effective modification of the Newton scheme consists in approximating the Jacobian by neglecting the contribution of the nonsymmetric matrix T , i.e. by setting: J ( E +G). This yields the Picard scheme, also known as successive iteration or nonlinear Richardson method [8]. There are a number of advantages of the Picard approach with respect to the Newton technique. First, the neglection of T leads to a symmetric system, with obvious efficiency gains in the linear solution phase. Second, the nonlinear behavior of Picard shows good initial convergence properties, especially in conjunction with relaxation techniques. On the other hand, the Picard approach often suffers from slow convergence or stagnation at relatively low residual levels. This stagnation can be put in connection with the neglection of the derivative terms of the Jacobian matrix. The initial solution estimate has a big effect on the behavior of iterative schemes. Unlike iterative techniques for linear systems, where theorems exist that guarantee global convergence under specified conditions, for Newton method only local convergence results are available, that is, provided the initial estimate is “close enough” to the solution. In practice, relaxed iteration, defined as um+1 = um + !sm (5) with ! 2]0; 1], can often help achieve or accelerate convergence when a poor initial estimate is used. 3.1 Solution of the Newton linearized system Newton equation (4) can be written as: 0B@ A Ad B C BT D0 +Dd 0 CT 0 0 1CA0B@ smq sm sm 1CA = 1 0BB@ Aqm; + B m; C m; + g BTqm; + f D0 m k CT m; 1CCA (6) where smq ; sm ; sm are the Newton directions for q; ; , respectively, D0 is defined as in the Picard iteration and Ad = @A @ qm; Dd = @D0 @ m; k Observing that matrices A, Ad, D0, Dd, and B are block-diagonal, the system (6) can be reduced to a system in the sm variable only. By eliminating smq and sm , the following system of equation is obtained: CTMCsm = CT 1 Mh+ ZH 1 hf D0( m k)i (7) where S = A 1B Z = A 1(B Ad) H = D0 +Dd +BTZ M = A 1 ZH 1ST and h = B m; C m; + g Nonsimmetric system (7) can be solved by preconditioned Krylov subspace methods [9]. The Newton directions for pressure head and velocity can then be obtained by sm = H 1 SCsm + Sh+ f D0( m k) (8) smq = 1 qm; + A 1 (B Ad)sm Csm + h (9) 4 Numerical Examples The first test case is a one-dimensional example in which the effect of an increasing gravity drainage component on Picard and Newton convergence is examined. The gravity drainage zone is defined as the region in a soil column which has zero pressure head gradient (and thus moisture flow is purely gravitational). As the length Lz of the soil column is progressively increased from 3 to 30 m, the extent of the gravity drainage component in the steady state infiltration solution also increases, as is seen in Figure 1 (left). The results of the simulation, summarized in Figure 1 (right), show that both Picard and Newton methods achieve convergence only if relaxed with very low values of ! (< 0:2). The mixed Picard-Newton scheme, where a number of initial relaxed Picard iterations are followed by full Newton iterations, converges rapidly. This behavior is explained by the fact that the initial Picard iterations provide an accurate initial solution for the subsequent Newton steps, which 0 10 20 30 40 50 60 70 80 90 100 # iterations 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 10 2 10 4 no nl in ea r er ro r no rm mixed (0.05) Picard (0.01) Newton (0.2) Newton (0.18) Newton (0.1−>1) Figure 1: Steady state solution and convergence profiles for test case 1. Scheme CPU # t Total # t time (s) nonl. it. min. max. avg. Picard 29.2 43 1000 4.8 10 3 1.00 0.12 Newton 5.3 28 101 6.0 10 3 0.91 0.18 Table 1: Convergence behavior of Picard and Newton scheme in the transient simulation. Test case 2a immediately achieve quadratic convergence. The Newton method with a relaxation parameter that increases towards unity with the progression of the iteration also achieves quadratic convergence, showing that line-search techniques can be effective in this example. The second test case involves two-dimensional transient flow in an unsaturated soil slab. The characteristic equations express the water saturation Sw in terms of effective saturation Se, in the form Sw( ) = (1 Swr)Se( ) + Swr, where Swr is the residual water saturation. The characteristic relations are then written as Se( ) = 8<: h1 + ( a ) i < a 1 a Kr( ) = Kr (Se( )) = Sen where a is the air entry pressure and , , , and n are constants. The parameter values used for the simulations are: = 0:5, x = 1 cm, z = 1 cm, Tmax = 5 d, (x; z; 0) = 60 cm. For test case 2a = 1, = 1, n = 1. For test case 2b = 3, = 3, n = 4. The characteristic feature of this example is the occurrence of convergence difficulties several time steps into the simulation, in response to the buildup of a sharp moisture front near the inflow boundary at x = 0, 6 z 10. Figure 2: Pressure head contours at time 5 days for test case 2b At 5 days the solution (shown in Figure 2) has essentially reached steady state. The convergence behavior of the Picard and Newton schemes is shown in Table 1. Note that at each time step of a transient simulation the nonlinear solver uses as initial guess the solution obtained at the previous time step. The difference between the initial and final solution of the nonlinear iteration is in general small. For this reason the Newton approach often achieves quadratic convergence at a very early stage of the iteration process. On the other hand, the Picard scheme still suffers from stagnation at a number of time steps. Convergence is achieved for all time steps only if relaxation with ! 0:4 is used, at the price of reduced efficiency. The steady state solution can also be obtained directly by solving the steady state equations. The convergence profiles of the nonlinear schemes are shown in Figure 3. For case 2a, the Picard scheme converges slowly with and without relaxation. The Newton and the mixed Picard-Newton methods share the same convergence profiles, with the mixed approach being slightly more computationally efficient. For case 2b, Picard converges only if relaxed (! = 0:4), and the mixed Picard-Newton technique achieves quadratic convergence after a halved number of iterations with respect to the full Newton approach, thus showing again the effectiveness of the mixed approach. 0 10 20 30 40 50 60 # iterations 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 10 2 10 4 no nl in ea r er ro r no rm Picard (0.4) Picard Newton mixed 0 10 20 30 40 50 60 # iterations 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 10 2 10 4 no nl in ea r er ro r no rm Picard (0.4) Picard Newton mixed Figure 3: Nonlinear convergence profiles for steady state simulation. Test case 2a (left) and 2b (right).